function [B,bf] = dirichlet_bc_bb(dofs,dim,V,E,T,ET,TE,cr,d,bdr,g,varargin)
% function [B,bf] = dirichlet_bc(dofs,dim,V,E,T,ET,TE,cr,d,bdr,g,varargin)
B = sparse(dim,dim); bf = zeros(dim,1);
nbdr = length(bdr);
J = (0:d)'/d; I = 1-J;
mat = vdm11(d);
for l = 1:nbdr
    eg = bdr(l); tri = ET(eg,1);
    eg2tri = find(TE(tri,:)==eg);
    v2_idx = mod(eg2tri,3)+1; 
    v3_idx = mod(mod(eg2tri,3)+1,3)+1;
    v2 = T(tri,v2_idx);
    v3 = T(tri,v3_idx);
    if eg2tri == 2
        eg2tri = -1* eg2tri;
    end
    % then i can get the edge's global dof;
    eg_dof = dofs(tri,cr_indices(0,d,eg2tri,cr)); 
    b_point_x = V(v2,1)*I + V(v3,1)*J;
    b_point_y = V(v2,2)*I + V(v3,2)*J;
%      plot(b_point_x,b_point_y,'r.');axis([-1 1 -1 1]); pause; 
    bf(eg_dof) = bf(eg_dof) + feval(g,b_point_x,b_point_y,varargin{:}); 
    B(eg_dof,eg_dof) = B(eg_dof,eg_dof) + mat;
end